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Nonlinear regression is a useful statistical tool, relating observed 
data and a nonlinear function of unknown parameters. When the 
parameter-dependent nonlinear function is computationally inten- 
sive, a straightforward regression analysis by maximum likelihood 
is not feasible. The method presented in this paper proposes to con- 
struct a faster running surrogate for such a computationally inten- 
sive nonlinear function, and to use it in a related nonlinear statis- 
tical model that accounts for the uncertainty associated with this 
surrogate. A pivotal quantity in the Earth's climate system is the 
climate sensitivity: the change in global temperature due to doubling 
of atmospheric CO2 concentrations. This, along with other climate 
parameters, are estimated by applying the statistical method devel- 
oped in this paper, where the computationally intensive nonlinear 
function is the MIT 2D climate model. 

1. Introduction. A fundamental question in understanding the Earth's 
climate system is quantifying the warming of the atmosphere due to in- 
creased greenhouse gases. This relationship is formalized by the climate 
sensitivity, a parameter defined as the increase in global mean surface tem- 
perature due to a doubling of CO2 in the atmosphere. Although climate 
sensitivity and other climate parameters are informed by observations, their 
impact can only be evaluated by simulations of climate with a numerical 
computer model. Such a model usually includes atmosphere, ocean, land and 
ice components and is called an atmosphere ocean general circulation model 
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(AOGCM). Because climate is denned as a long term average of weather, an 
AOGCM is usually run (integrated) over many years in order to establish 
its mean behavior. Thus, numerical experiments with these models require 
extensive computational resources and the number of runs (also termed in- 
tegrations) is often limited. For example, the Community Climate System 
Model (CCSM) requires months of time on a supercomputer to simulate 
a few hundred model years. Typically an AOGCM depends on unknown 
parameters which need to be estimated and the statistical challenge is to 
estimate these parameters along with companion measures of uncertainty 
using a limited set of climate model experiments. 

The general statistical problem addressed here is parameter estimation 
in a nonlinear regression model [e.g., Seber and Wild (1989)], where the 
nonlinear regression function is computationally intensive to evaluate, such 
as an AOGCM. In particular, the example discussed here involves observed 
climate data and the MIT 2D climate model [Sokolov and Stone (1988)] 
as a nonlinear function of three uncertain parameters collectively denoted 
by 9: the equilibrium climate sensitivity S, the diffusion of heat anomalies 
into the deep ocean K v and the net aerosol forcing F aer . Here we take a 
maximum likelihood approach and pay particular attention to the correla- 
tion structure and its effects on the uncertainty measures of the resulting 
estimates. The standard nonlinear regression approach for this particular 
estimation problem assumes the following statistical model for the observed 
climate data: 

(1) Y = f e + e, 

where the errors are assumed normal with zero mean vector and covariance 
matrix W. The estimated parameters of this statistical model, including 9, 
are then obtained by maximizing the likelihood 

< 2) (ir^ exp B (Y - w,F ~ i(Y -4 

This is usually achieved by using an iterative algorithm. Notice, however, 
that this requires computing fg for many values of 9, or, equivalently, run- 
ning the climate model for a possibly large number of 9 values. Such an 
approach is not feasible for applications where fg is an AOGCM or even the 
simplified MIT 2D climate model. To overcome this computational difficulty, 
we will substitute a statistical surrogate for fg, denoted fg, which will result 
in a much faster estimation algorithm for the unknown parameters 9. The 
statistical model for the observed data is now 

(3) Y = ~fg + B e + e, 

where E# is the error in approximating fg by fg and e is observation error. 
More precisely, we use analysis of computer experiments methodology [e.g., 



PARAMETER ESTIMATION FOR NONLINEAR REGRESSION 



3 



Santner et al. (2003)] to analyze climate model output data at a sample 
of 'input' parameter vectors and build a statistical surrogate to predict the 
climate model output at new, untried parameter values. This analysis is 
empirical Bayesian in its nature [as described in Currin et al. (1991) for 
scalar output and in Drignei (2006) for multiple outputs], in the sense that 
a Gaussian process serves as a prior distribution for the multiple outputs and 
the posterior distribution is used to predict the output at new parameters. 
We take to be the posterior mean, but also use the uncertainty about 
this mean to adjust the likelihood of observations. To our knowledge, the 
empirical Bayesian analysis of multiple output computer experiments, used 
in the context of physical model calibration through maximum likelihood 
methods, is new. The statistical model we develop also accounts for different 
types of uncertainty (e.g., climate model internal variability), it includes 
correlation with respect to various dimensions (e.g., space-time correlation) 
and a new feature of our work is to include the uncertainty of the surrogate 
model as part of the estimated parameter uncertainty. The calibration of 
computationally intensive computer models has been recently done mostly 
through fully Bayesian models, for example, Kennedy and O'Hagan (2001), 
Craig et al. (2001), Bayarri et al. (2007), Higdon et al. (2008) and the last 
two references dealing with multidimensional computer model outputs. 

We make use of the data sets from Forest et al. (2002) who implemented a 
Bayesian model for the same problem and have shown results from flat and 
expert priors (see top row of Figure 2). A notable feature from a climatolog- 
ical point of view in their work is the skewness of the climate sensitivity S, 
which appears to be more pronounced for flat priors. Appendix A presents a 
geophysical argument in support of this skewness. Forest et al. (2003), For- 
est, Stone and Sokolov (2006) and Sanso et al. (2008) revisited and refined 
the Bayesian approach. Several other studies estimated a probability density 
function for climate sensitivity [Andronova and Schlesinger (2001), Gregory 
et al. (2002), Knutti et al. (2002, 2003)]. All such studies are based on es- 
timating the degree to which a climate model can reproduce the historical 
climate record. More recently, this same technique has also been applied to 
the climate record for the past 600 years [Hegerl et al. (2006)]. Last Glacial 
Maximum (LGM) climate change can also be used as an additional line of 
evidence to estimate the probability density function of the climate sensitiv- 
ity and such results have been combined by Annan and Hargreaves (2006) to 
provide a more complete picture of available constraints for placing bounds 
on climate sensitivity. 

This paper is organized as follows. Section 2 discusses the climate ob- 
servations, the climate model and the output data sets. Section 3 develops 
the statistical surrogate for the climate model, while Section 4 shows how 
this surrogate can be used to construct a computationally efficient statistical 
model for the climate observations. The results are presented in Section 5 
and the paper ends with some conclusions in Section 6. 
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2. Observed and modeled climate. 

2.1. Climate observations. In this analysis we use three separate sets of 
observations representing the ocean heat content, surface temperatures and 
upper air temperatures. 

The oceans play an important role in the planet's climate system because 
they can store and transport large amounts of heat. The subsurface ocean 
temperature records are sparse and contain the most uncertainty due pri- 
marily to the difficulty in obtaining such temperature data sets. The data 
used in this paper originate in the Levitus et al. (2000) ocean temperature 
data set from the surface through 3000-meter depth, showing a net warming 
over the last half of the twentieth century. 

The upper-air temperature data set used in this paper originates in the 
improved upper-air temperature data recorded by radiosondes and described 
in Parker et al. (1997). The later data are considered somewhat more re- 
liable than satellite-based Microwave Sounding Units (MSU) data because 
they provide a longer record and a better vertical resolution. Nevertheless, 
the MSU data have proved to be useful in removing time-varying biases of 
radiosonde temperature data caused, for example, by changes of instrumen- 
tation or operating procedures. Parker et al. (1997) describe some interpola- 
tion, bias- and error-removing methods for the original sparse and irregular 
radiosonde upper-air temperature data sets. 

Among all the temperature records, the surface temperatures are the 
longest, most spatially complete and documented. This is due primarily to 
the existence of meteorological stations throughout the world for a relatively 
long time. The data sets used in this paper originate in the extended and 
interpolated data sets of Jones (1994). 

Summary statistics for each observational dataset are used to make the 
appropriate comparison with the climate model (as discussed later). The 
averaging methods for each diagnostic are discussed in Forest et al. (2002) 
and serve to remove the short time-scale variability that is not associated 
with the long time-scales changes in climate. These averages result in the 
patterns summarizing the changes in mean temperatures for each source as 
discussed in Section 2.3. 

2.2. The MIT 2D climate model and the unknown parameters. AOGCMs 
are the primary tools for predicting changes in global climate patterns on 
planetary scales and at decadal or longer time scales. Mathematically, these 
models are systems of partial differential equations derived from the laws of 
fluid dynamics and thermodynamics for the atmosphere, ocean, ice and land 
systems, in three spatial dimensions (3D) and a temporal dimension. Since 
these models are nonlinear, they are usually solved numerically over a space- 
time grid. These numerical methods, however, require large computational 
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resources and, therefore, have limited flexibility for exploring parametric 
(or structural) uncertainty. Most often, for a given AOGCM, individual pa- 
rameters are set according to performance of individual components (e.g., 
cloud parameterizations or ocean sub-grid scale mixing parameterizations) 
and then held fixed or modified in a heuristic fashion when the fully coupled 
AOGCM is assembled and tested. For a more statistical approach to deter- 
mine parameters, we require a flexible climate model designed to explore 
uncertainty in the large-scale response (e.g., global or hemispheric average 
temperature changes) in the fully coupled system. The MIT 2D climate 
model [Sokolov and Stone (1988)] is one such model designed for these pur- 
poses and was used in this research. Some technical aspects of the MIT 2D 
climate model are given in Appendix B. 

The MIT 2D model has two parameters that determine the decadal to 
century response to external factors that drive the climate. These are the 
equilibrium climate sensitivity S (measured in degrees K) to a doubling 
of CO2 concentrations and the global-mean vertical thermal diffusivity K v 
(measured in cm 2 /s) for the mixing of thermal anomalies into the deep 
ocean. Sokolov and Stone (1988) have shown that the large-scale response 
of AOGCMs can be duplicated by the MIT 2D model with an appropriate 
choice of these two parameters. This correspondence supports the study of 
the climate system using the simpler MIT model because it can approximate 
more complex and realistic three dimensional models. The third parameter 
considered in this paper is the net aerosol forcing strength F aeT (measured in 
W/m 2 ) and controls the amount of cooling of the atmosphere to increased 
amounts of particles. The computational time required to simulate 50 years 
with the MIT 2D climate model is about 4 hours on a 3 GHz Pentium 4 Linux 
workstation and it is several orders of magnitude faster than simulation with 
a state-of-the-art 3D AOGCM. Given its flexibility to duplicate AOGCM 
responses and its computational efficiency, the MIT 2D climate model is 
a tool well suited for answering questions which would be impractical to 
explore with 3D AOGCMs. 

2.3. Specific data sets and model output. Due to the sparsity and the 
large uncertainties in the deep-ocean temperature data, only the linear trend 
in the observed temperature record is retained for analysis so that the ocean 
observed data is just a scalar. The upper-air temperature observations are 
the difference in the 1986-1995 and 1961-1980 mean temperatures, recorded 
at each 5° latitude and at 8 pressure levels (850 hPa through 50 hPa). In 
order to simplify the analysis, 10 latitude coordinates containing mostly 
missing data have been discarded, so that the final upper-air temperature 
change observed data set is a 26 x 8 matrix. The surface temperature data 
set has been averaged so that the final surface temperature change data set 
is a 4 x 5 matrix, corresponding to 4 latitude bands by 5 decades. 
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FlG. 1. The sampled inputs. 



Let 6 = [S, Kj , F acr ] be the parameter vector. We considered D = 306 
parameters 6^ as in Forest et al. (2002), sampled in a parameter space (see 
Figure 1) to cover a range of parameters well beyond the domain correspond- 
ing to current AOGCMs. In Forest et al. (2002), an initial, approximately 
factorial design over the domain has been selected and then a second, higher 
density sampling has been added in the region of highest likelihood to better 
estimate their entire distribution. For each parameter vector in this param- 
eter space, the relatively large model output data set is transformed so that 
its format matches that of the observed data sets: the deep-ocean temper- 
ature trend (scalar), the upper-air temperature changes (26 latitudes and 8 
levels) and the surface temperature changes (4 latitude zones and the last 5 
decades of the last century) . There are four replicates for each of these data 
sets, called ensemble members, obtained by changing the initial conditions 
in the climate model. In this paper we work with averages across ensem- 
ble members, although the ensemble variability will be accounted for in our 
statistical model. 

3. Statistical surrogate for the climate model. The statistical model for 
observations developed in this paper has two components. The first compo- 
nent is a surrogate of the climate model output over the region of model 
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parameters and the second component connects the surrogate with the ob- 
served data. Each of these components will be discussed next. 

3.1. Statistical model for climate model output. All three observational 
data sets have a similar statistical model, so we will only present the model 
for surface temperature changes in detail. The climate model surface tem- 
perature change can be conveniently organized as an array over time, space 
and climate parameters, with dimension Nt x Nz x D. (For our data sets, 
Nt = 5, Nz = 4, D = 306.) To better describe the statistical model, this 
three-dimensional data set is stacked as a vector and simply denoted by f . 
Climatology arguments suggest that the numerical climate model be decom- 
posed into two random components: f = f s + f n corresponding to a climate 
signal f s and noise process f n ascribed to climate model internal variability. 

The correlation matrix for the signal f s can be written, under the compu- 
tationally convenient assumption of separability, as a Kronecker product of 
smaller correlation matrices Cq <g> Cz <8> Ct reflecting the various dimensions 
of the problem. The power exponential [Sacks et al. (1989)] was chosen as a 
simple but flexible model to describe the correlations among each dimension 
of climate signal. The matrix Ce is defined as 



The correlations among zonal bands and across time are also assumed to 
follow the power exponential family, although analysis based on the Matern 
family [e.g., Stein (1999), page 31] has also been carried out. The correlation 
matrix for the noise f n will be I <g> T, with a consistent empirical estimator 
of the space—time covariance T obtained by assuming the same space-time 
structure across the D sampled parameters and 4 ensemble members (after 
subtracting the ensemble mean at each of the D sampled parameters). A 
simple method based on wavelet basis functions [Nychka et al. (2002)] was 
used to estimate a nonstationary form that takes advantage of the assump- 
tion that r does not depend on 9 and can be thought of as a regularization 
of the sample covariance matrix. The space time fields are expanded in a W- 
transform multiresolution basis, with W being the matrix of wavelet basis 
functions that relate coefficients to the model field and T> being the sam- 
ple covariance matrix for the coefficients. T> is decomposed as T> = TC 2 and 
the elements of 7i are thresholded by setting 90% of the elements to zero, 
resulting in a regularized matrix Ti.. We use the estimate T = W7i. 2 W . Then 



is the overall covariance matrix, so that f ~N(/il,Se). The output data 
analyzed here are temperature changes and averages, therefore, we assume 





S© = a 2 (C e ®C Z ® C T ) + u 2 (I <g> T) 



<s 
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that any large scale trends have been eliminated, justifying the choice of a 
statistical model with constant mean. An unbiased estimate of the parame- 
ter uj 2 is obtained by pooling the ensemble sample variances across inputs, 
space and time. The remaining statistical parameters are estimated by the 
maximum likelihood method and all parameters are fixed at their estimated 
values throughout the rest of the statistical analysis. 

3.2. The statistical surrogate and its error. The statistical model de- 
scribed above is used to construct a surrogate for the climate model at an 
arbitrary parameter vector in the parameter space. To define the surrogate 
for fg and its error, we consider the conditional distribution of the climate 
signal f| on the climate model output data f , which is multivariate normal 
of mean vector 

fa = fj,l + See^eHf ~ A* 1 ) 

and covariance matrix 

V e = a 2 (C z ® C T ) + co 2 T - E^Eg 1 ^ 

where = ct 2 (Cqq <g> Cz <8> Ct) and Cqq is defined as in (4). The surrogate 
is {q and its error E has a multivariate normal distribution N(0, Vq). The 
surrogate model for the upper air temperature is defined similarly. The sur- 
rogate model for the deep ocean temperature trend is much simpler because 
it is based on a univariate response. 

4. Likelihood function for the observed data. 

4.1. Likelihood function for a single data set. Here we develop the model 
for the observed surface temperatures, whereas Section 4.2 presents the com- 
plete statistical model for all three data sets. 

The statistical surrogate developed above will not predict perfectly the 
climate model output at new parameters, nor do we expect it to match 
perfectly the observed data. Therefore, if Y denotes generically the observed 
data set in vector format, then the corresponding statistical model is 

Y = f e + E e + e, 

where e is observation error with multivariate normal distribution N(0, t 2 Rz® 
Rt) and the elements of the matrices Rz and Rt are power exponential cor- 
relations with parameter vector Here E and e are assumed independent. 
An important statistical point is that the discrepancy between f# and fg is 
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modeled explicitly through the analysis of the computer experiments' ap- 
proach in the previous section. The likelihood of the model for Y is 



Whereas the computationally intractable likelihood (2) contains the non- 
linear climate model fg and the covariance matrix W of the observation error, 
the more computationally efficient likelihood (5) depends on the surrogate 
fg of the climate model and a covariance matrix based on two components: 
surrogate error and observation error. 

4.2. Full likelihood. Let Y S ,Y ,Y U denote the observed surface temper- 
ature change, the observed deep ocean temperature trend and the observed 
upper air temperature change respectively. These three data sets are as- 
sumed independent conditioned on the true climate signal at 9 so that the 
overall likelihood to be optimized is a product of the likelihood functions 
developed in the previous subsection, for the three specific data sets: 



L(Y s ,Y ,Y u \e,T s ,T ,T u ^ s ,Cu)=L(Y s \6,T s ,Cs)L(Y \9,T )L(Y u \9,T u ,Cu)- 



(Notice that the observation error for the ocean temperature trend is uni- 
variate normal.) This likelihood takes about 5 seconds to be evaluated in 
Matlab on a computer with dual 2.6 GHz Xeon processors and 4 GB RAM. 

There are geophysics arguments in support of a right skewed distribution 
for the climate sensitivity parameter (see Appendix A) and the literature on 
the estimation of this parameter reports various degrees of skewness. This, in 
turn, makes it difficult to agree on an upper confidence bound for this impor- 
tant parameter: how large will the global mean surface temperature be when 
the CO2 amount is doubled? Here we would like to investigate this aspect 
through the finite sample distribution of the maximum likelihood estimators. 
An established method for such purpose is the parametric bootstrap [Efron 
and Tibshirani (1993)], where synthetic data Y S ,Y D ,Y U will be generated 

from the multivariate normal of likelihood L(Y s ,Y ,Y u \9,f s ,f ,f u ,^ s ,^ u ), 
with the maximum likelihood estimates taken as "true" values of the param- 
eters. The new likelihood of the simulated data will be maximized and the 
point estimate (#,r,£) will be obtained. This process is repeated B times 
and, therefore, B independent simulated values (9, f,£) will be obtained, 
which will further be used to summarize the distribution of (#,f,£), for ex- 
ample, through confidence regions. A direct search method in Matlab has 
been used to optimize the likelihood functions throughout this paper. 



(5) 





10 



D. DRIGNEI, C. E. FOREST AND D. NYCHKA 



5. Results. The statistical model described in Section 3.1 has been fitted 
to the output data and the parameter estimates for the power exponential 
correlations given in Table 1. An important component of the statistical 
model is the choice of covariance family for the statistical surrogate. To 
determine the sensitivity to the power exponential correlation, the Matern 
family was also considered and we present how the inference on the climate 
parameters is effected by this alternative model for correlations. Table 1 



Table 1 

Estimates of parameters in the statistical model for output data (the power exponential 
model): surface temperature (upper row), upper air (middle row) and deep ocean (lower 

row) 



M f?i Pi r\i Pi f?3 P3 




Pt Vs 


Ps 


a 2 


u, 2 


0.217 1.365 0.425 1.189 0.273 2.283 0.903 
-0.047 1.206 0.302 1.031 0.185 1.274 0.390 
0.001 3.430 1.999 21.636 1.999 1.545 1.927 


1.227 
14.476 


1.147 1.255 
1.849 4.382 


1.501 
1.431 


0.024 
0.007 

2 x 10~ 6 


0.016 
0.004 

7 x 10~ 9 




02468 10 2 4^ a -1.5 -1 -0.5 



Fig. 2. Upper row: estimated univariate marginal densities based on power exponential 
correlations (solid) with 99% confidence intervals (diamond), and on the Matern correla- 
tions (dot) with 99% confidence intervals (square), and the posterior pdfs from Forest et al. 
(2002): expert prior for climate sensitivity (dash) and flat priors (dash-dot); Lower row: 
estimated bivariate densities with contours representing 99%, 95% and 90% confidence 
regions (solid lines: power exponential; dotted lines: Matern). 
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Table 2 

Estimates of parameters in the statistical model for output data (the Matern model): 
surface temperature (upper row), upper air (lower row). The deep ocean estimates and u> 2 

are the same as in Table 1 





Vi 


Pi 


V2 


P2 


V3 


P3 


a t 


a s 


<r 2 


0.211 
-0.021 


1.399 
1.142 


0.426 
0.323 


1.239 
0.963 


0.267 
0.193 


2.287 
1.171 


0.908 
0.376 


3.987 
6.873 


2.682 
12.244 


0.022 
0.006 



indicates smoothness in the latitude-time dimensions for surface tempera- 
tures and latitude-pressure dimensions for upper-air temperatures. There- 
fore, Matern correlations generating one-time mean square differentiable re- 
alizations (i.e., with the smoothness parameter v set to 1.5) have also been 
fitted in those dimensions, with the correlation parameter estimates shown 
in Table 2. More specifically, the statistical model for the output data based 
on the Matern correlations is identical to the model presented in Section 3, 
except that the elements of the matrices Cz and Ct are now Matern instead 
of power exponential correlations, with temporal and latitude range param- 
eters denoted by at and a s , respectively. The power exponential maximized 
loglikelihood minus the Matern maximized loglikelihood is 97 for the surface 
temperature and 1555 for the upper air output data sets, respectively. The 
models with these correlations have been further used to build the surro- 
gates and the associated errors as described in Section 3.2. These have been 
used in the likelihood of the observed data (Section 4). B = 300 bootstrap 
simulations were obtained, which have been further used to obtain nonpara- 
metric kernel marginal density estimates and confidence intervals/regions. 
The results are summarized in Figure 2. The upper row contains the univari- 
ate marginal density estimates with 99% confidence intervals for both power 
exponential (solid) and Matern (dot) models, as well as posterior pdfs from 
Forest et al. (2002); it appears that a mild right skewness of the estimated 
densities for the climate sensitivity S is present. The lower row contains the 
kernel density estimates of the bivariate marginal densities under the power 
exponential (solid) and Matern (dot) models. The contours correspond to 
90%, 95% and 99% levels confidence regions. The results from these two cor- 
relation models (power exponential and Matern) do not appear to be largely 
different. The inclusion of correlations was helpful in better constraining the 
climate system parameters than in the statistical model considered in Forest 
et al. (2002), which did not include such correlations. 

6. Conclusion. This paper has presented a computationally efficient sta- 
tistical model as an alternative to a naive nonlinear regression model in- 
volving a computationally intensive nonlinear function. Our approach is to 
construct a statistical model that includes a faster running surrogate for the 
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computationally intensive nonlinear function and it is underlined by some 
general principles. The statistical model offers a comprehensive framework to 
accommodate various sources of uncertainty. For example, we incorporated 
an error term to represent climate internal variability and a component for 
uncertainty in the surrogate, as well as an observation error term. The co- 
variance matrices are full rank and estimated from the available forced model 
output data. The statistical model proposed accounts for correlation in all 
dimensions of the problem (e.g., it includes spatio-temporal correlation). 
This strategy was helpful to better constrain the unknown climate model 
parameters through tighter confidence intervals and regions, especially for 
the climate sensitivity S and ocean diffusivity K v parameters. 

APPENDIX A: CLIMATE SENSITIVITY 

An important characteristic in the Earth's climate system is the climate 
sensitivity S: the change in global temperature due to doubling of atmo- 
spheric CO2 concentrations. In complex climate models derived from the 
fundamental equations of physics, the individual processes interact and the 
feedbacks between the key processes lead to the overall amplification of the 
temperature change. At the basic level, the additional greenhouse gases (e.g., 
carbon dioxide, methane and nitrous oxide) have a direct impact on the out- 
going longwave radiation in the atmosphere and in the absence of feedbacks 
(i.e., the nontemperature atmospheric state variables remains unchanged), 
this change to the system will eventually lead to a given amount of tem- 
perature change. But, when the temperature changes, this will inevitably 
lead to additional changes in the temperature and humidity profiles, cloud 
distributions, snow and ice cover, and other state variables. These additional 
changes are the climate system feedbacks that can amplify or reduce the di- 
rect impact of the greenhouse gases on the radiative forcing. If we let AT Q be 
the temperature change due to the direct impact on the radiative transfer, 
we note that AT D includes one feedback, the direct impact of temperature 
change on the infrared spectrum by changing greenhouse gas concentrations, 
but does not include additional feedbacks from other components of the cli- 
mate system. It is customary to let the total change in global mean surface 
temperature above the equilibrium value be written as AT2 X = [e.g., 
see Hansen et al. (1984) or Schlesinger and Mitchell (1987)]. Thus, we have a 
one-to-one correspondence between the feedback and the climate sensitivity 
(AT2 X and the uncertainties in each can be related). 

This relation between climate sensitivity and feedbacks has been known 
since the early days of climate research and was described in the Charney 
Report [NRC (1979)], formalized in Hansen et al. (1984) and included in 
textbooks as early as Henderson-Sellers and McGuffie (1987). Hence, we 
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now have a simple expression to relate the uncertainty in these three vari- 
ables and from first principles, AT Q and 4> are the two variables which can be 
quantified individually using climate models and combined to yield AT2 X - 
Schlesinger and Mitchell (1987) quantified the uncertainties in feedbacks and 
recognized that the uncertainties in individual processes will combine to pro- 
vide an estimate of the total uncertainty. Given the number of processes in 
the system, it is typical to consider the net feedback to have a normal dis- 
tribution. However, we also need to consider AT Q as an additional uncertain 
variable with its own normal distribution, although this is a small contri- 
bution compared to (j). By combining these two distributions for AT Q and 
(j), the resulting distribution for climate sensitivity is expected to be right 
skewed, although one cannot clearly infer only from these arguments how 
long the right tail will be. 

APPENDIX B: TECHNICAL ASPECTS OF THE MIT 2D CLIMATE 

MODEL 

The MIT 2D climate model consists of a zonally averaged atmospheric 
model coupled to a mixed-layer Q-flux ocean model, with heat anomalies 
diffused below the mixed-layer. The model details can be found in Sokolov 
and Stone (1988). The atmospheric model is a zonally averaged version of 
the Goddard Institute for Space Studies (GISS) Model II general circulation 
model [Hansen et al. (1983)] with parameterizations of the eddy transports 
of momentum, heat and moisture by baroclinic eddies [Stone and Yao (1987, 
1990)]. The model version we use has 24 latitude bands and 11 vertical layers 
with 4 layers above the tropopause. 

The model also employs a Q-flux ocean mixed layer model with diffusion 
of heat anomalies into the deep-ocean below the climatological mixed layer. 
This model of the ocean component of the climate system is fully described 
by Hansen et al. (1983) and only increased computations by a few percent. 
The model uses the GISS radiative transfer code which contains all radia- 
tively important trace g well as aerosols and their effect on radiative 

transfer. The surface area of each latitude band is divided into a percent- 
age of land, ocean, land-ice, and sea-ice with the surface fluxes computed 
separately for each surface type. This allows for appropriate treatment of ra- 
diative forcings dependent on underlying surface type such as anthropogenic 
aerosols. The atmospheric component of the model, therefore, provides most 
important nonlinear interactions between components of the atmospheric 
system. 
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